﻿#include "raster_surface.h"

#include <fstream>
#include <algorithm>
using namespace std;
namespace oil
{
RasterSurface::RasterSurface(const Vector4f& min, const Vector4f& max,
	float width, float height) :min_(min), cell_width_(width), cell_height_(height)
{
	int rows, cols;
	rows = std::ceil((max.y() - min.y()) / cell_height_) + 1;
	cols = std::ceil((max.x() - min.x()) / cell_width_) + 1;

	rows_ = rows;
	cols_ = cols;

	data_.resize(rows * cols);

}

RasterSurface::RasterSurface(const Vector4f& min, int rows, int cols,
	float width, float height) :min_(min), rows_(rows), cols_(cols),
	cell_width_(width), cell_height_(height)
{
	data_.resize(rows * cols);
}

RasterSurface::RasterSurface()
{
	rows_ = 0;
	cols_ = 0;
	min_ = Vector4f::Zero();
	cell_width_ = 0;
	cell_height_ = 0;
	
}

RasterSurface::~RasterSurface()
{
}
double RasterSurface::operator()(int row, int col) const
{
	int pos_r = row < 0 ? 0 : row;
	pos_r = pos_r > rows_ - 1 ? rows_ - 1 : pos_r;
	int pos_c = col < 0 ? 0 : col;
	pos_c = pos_c > cols_ - 1 ? cols_ - 1 : pos_c;
	int pos = pos_r * cols_ + pos_c;
	return data_[pos];
}

double RasterSurface::operator()(double x, double y)
{
	int row = std::floor((y - min_.y()) / cell_height_);
	int col = std::floor((x - min_.x()) / cell_width_);

	//获取4X4 格网的数据
	double x_0_1 = std::fmod(x - min_.x(), (double)cell_width_) / cell_width_;
	double y_0_1 = std::fmod(y - min_.y(), (double)cell_height_) / cell_height_;
	double p[4][4];
	for (int r = -1; r < 3; ++r)
	{
		for (int c = -1; c < 3; ++c)
		{
			int pos_r = row + r;
			int pos_c = col + c;
			pos_r = pos_r < 0 ? 0 : pos_r;
			pos_r = pos_r > rows_ - 1 ? rows_ - 1 : pos_r;
			pos_c = pos_c < 0 ? 0 : pos_c;
			pos_c = pos_c > cols_ - 1 ? cols_ - 1 : pos_c;

			int pos = pos_r * cols_ + pos_c;

			p[r + 1][c + 1] = data_[pos];
		}
	}
	double z = bicubicInterpolate(p, x_0_1, y_0_1);

	return z;
}

float& RasterSurface::operator()(int row, int col)
{
	int pos_r = row < 0 ? 0 : row;
	pos_r = pos_r > rows_ - 1 ? rows_ - 1 : pos_r;
	int pos_c = col < 0 ? 0 : col;
	pos_c = pos_c > cols_ - 1 ? cols_ - 1 : pos_c;
	int pos = pos_r * cols_ + pos_c;
	return data_[pos];
}

double RasterSurface::cubicInterpolate(double p[4], double x) const
{
	return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] +
		x*(3.0*(p[1] - p[2]) + p[3] - p[0])));
}

double RasterSurface::bicubicInterpolate(double p[4][4], double x, double y) const
{
	double arr[4];
	arr[0] = cubicInterpolate(p[0], y);
	arr[1] = cubicInterpolate(p[1], y);
	arr[2] = cubicInterpolate(p[2], y);
	arr[3] = cubicInterpolate(p[3], y);
	return cubicInterpolate(arr, x);
}

void RasterSurface::save(const string& file_path)
{
	ofstream fout(file_path);

	auto minmax = minmax_element(begin(data_), end(data_));
	float maxX = min_.x() + cols_ * cell_width_;
	float maxY = min_.y() + rows_ * cell_height_;

	if (fout.is_open())
	{
		fout << "DSAA" << endl;

		fout << cols_ << " " << rows_ << endl;
		fout << min_.x() << " " << maxX << endl;
		fout << min_.y() << " " << maxY << endl;
		fout << *minmax.first << " " << *minmax.second << endl;

		for (int i = 0; i < rows_; ++i)
		{
			for (int j = 0; j < cols_; ++j)
			{
				int pos = i * cols_ + j;
				fout << data_[pos] << " ";
			}
			fout << endl;
		}
	}
}
}//end oil
